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The energy of a free gas of neutrons and protons is well known to be approximately isospin 
parabolic with a negligibly small quartic term of only 0.45 MeV at the saturation density of nuclear 
matter po = 0.16 fm -3 . Using an isospin-dependent single-nucleon momentum distribution including 
a high (low) momentum tail (depletion) with its shape parameters constrained by recent high-energy 
electron scattering and medium-energy nuclear photodisintegration experiments as well as the state- 
of-the-art calculations of the deuteron wave function and the equation of state of pure neutron matter 
near the unitary limit within several modern microscopic many-body theories, we show for the first 
time that the kinetic energy of interacting nucleons in neutron-rich nucleonic matter has a significant 
quartic term of 7.18 ± 2.52 MeV. Such a large quartic term has broad ramifications in determining 
the equation of state of neutron-rich nucleonic matter using observables of nuclear reactions and 
neutron stars. 
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I. INTRODUCTION 

To determine the equation of state (EoS) of isospin- 
asymmetric nuclear matter (ANM) has been a long¬ 
standing goal shared by both nuclear physics and as¬ 
trophysics [1]. Usually one uses the so-called empirical 
parabolic law for the energy per nucleon, i.e., E(p,S) = 
E 0 (p) + E sym (p)S 2 + 0(S 4 ) where p = p n + p p and 
5 = {p n — Pp)/P are the nucleon density and isospin asym¬ 
metry of the system in terms of the neutron and proton 
densities p n and p p , respectively. The isospin quadratics 
of the ANM EoS has been verified to high accuracies from 
symmetric (5 = 0) up to pure neutron (d = 1) matter by 
most of the available nuclear many-body theories using 
various interactions, see, e.g., ref. [2], Nevertheless, it has 
been shown consistently in a number of studies that for 
some physical quantities relevant for understanding prop¬ 
erties of neutron stars, such as the proton fraction at /3 
equilibrium, core-crust transition density and the criti¬ 
cal density for the direct URCA process to happen, even 
a very small coefficient E sym ^(p) of the isospin quartic 
term in the EoS can make a big difference [3]. 

Here we concentrate on examining the isospin quadrat¬ 
ics of the kinetic EoS. For many purposes in both nu¬ 
clear physics and astrophysics, such as simulating heavy- 
ion collisions [4] and determining critical formation densi¬ 
ties of different charge states of A resonances in neutron 
stars [5], one has to know separately the kinetic and po¬ 
tential parts of the EoS. While neither any fundamental 
physical principle nor the empirical parabolic law of the 
EoS requires the kinetic and potential parts of the EoS 
to be quadratic in S individually, in practice especially in 
most phenomenological models the free Fermi gas (FFG) 
EoS is often used for the kinetic part and then the gen¬ 
erally less known potential EoS is explored by compar¬ 
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ing model predictions with experimental data. It is well 
known that the FFG model predicts a kinetic symme¬ 
try energy of if^(po) ~ 12.3 MeV and a negligibly 
small quartic term of E^ 4 (po) = E^ , f a (p 0 )/27 « 0.45 
MeV at po = 0.16/fm 3 . However, nuclear interactions, 
in particular the short-range repulsive core and tensor 
force, lead to a high (low) momentum tail (depletion) in 
the single-nucleon momentum distribution above (below) 
the nucleon Fermi surface [6-9]. Much progress has been 
made recently both theoretically and experimentally in 
quantifying especially the nucleon high momentum tails 
(HMT) in ANM, see, e.g., refs. [10-14], In this work, us¬ 
ing isospin-dependent nucleon HMT constrained by re¬ 
cent high-energy electron scattering and medium-energy 
nuclear photodisintegration experiments as well as the 
state-of-the-art calculations of the deuteron wave func¬ 
tion and the EoS of pure neutron matter (PNM) near the 
unitary limit within several modern microscopic many- 
body theories, we show that the kinetic ANM EoS has a 
significant quartic term of E^^^po) = 7.18 ± 2.52 MeV 
that is about 16 times the FFG model prediction. 

II. ISOSPIN DEPENDENCE OF 
SINGLE-NUCLEON MOMENTUM 
DISTRIBUTION WITH A HIGH MOMENTUM 
TAIL IN NEUTRON-RICH MATTER 

Guided by well-known predictions of microscopic nu¬ 
clear many-body theories, see, e.g., reviews in ref. [15], 
and recent experimental findings [10-13], we describe the 
single-nucleon momentum distribution in ANM using 

( A J+ pjl(\k\/kJ), 0 < |k| < kp, 

K{p^) = { 

l Cj (fcp/|k|) , fcp < |k| < </>jfcp. 

(1) 

Here, J = n,p is the isospin index, fcp = fcp(l + r/d) 1 ' 3 
is the transition momentum [13] where fcp = (3tt 2 p/2) 1 ^ 3 
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where M is the nucleon mass. To ensure that the 
i?j cm (p) vanishes as required by the neutron-proton ex¬ 
change symmetry of the EoS, we require that A" = — A^, 
/3" = — Pi, Cf = —Cf and 0" = —0f, i-e., more com¬ 
pactly Yj = F 0 (l + Air/ 5). 


FIG. 1: A sketch of the single-nucleon momentum distribution 
with a high momentum tail. 

are depicted in Fig. 1. The A j measures the depletion 
of the Fermi sphere at zero momentum with respect to 
the FFG model prediction while the /3j is the strength 
of the momentum dependence J(|k|/fc/) [16-18] of the 
depletion near the Fermi surface. The jump Z / of the 
momentum distribution at fc/, namely, the “renormal¬ 
ization function”, contains information about the nu¬ 
cleon effective E-mass and its isospin dependence[19]. 
Specifically, = n^j _ 0 - n/j +0 = where 

= [1 — dV /9w ] _1 with V and ui being the real 
part of the single-particle potential and energy [ 6 , 20 ], 
respectively. 

The amplitude Cj and cutoff coefficient (f>j determine 
the fraction of nucleons in the HMT via 

A MT = 3 cj (1 - A). (2) 

The normalization condition [2/ ( 27 t) 3 ] n£(p, d)dk = 

Pj = (fc/) 3 / 37 r 2 requires that only three of the four pa¬ 
rameters, i.e., /3 j, Cj, <f)j and Aj, are independent. Here 
we choose the first three as independent and determine 
the A j from 


Aj = 1 — 


3 Pj [ ki 

(^f ) 3 Jo 



k 2 dk - 3 Cj 



(3) 

Hinted by the finding within the self-consistent Green 
function (SCGF) theory [21] and the Brueckner-Hartree- 
Fock (BHF) theory [22] the depletion A j has an almost 
linear dependence on 5 in the opposite directions for neu¬ 
trons and protons, we expand all four parameters in the 
form Yj = Fo(l + A/d). Then, the total kinetic energy 
per nucleon in ANM 


E kin ( P ,6) 




(p, d)dk 


(4) 


III. CONSTRAINING THE PARAMETERS OF 
THE SINGLE-NUCLEON MOMENTUM 
DISTRIBUTION 

It is well known that the nucleon HMT from deuteron 
to infinite nuclear matter scales, see, e.g., refs. [23-25], 
leading to constant per nucleon inclusive (e, e') cross sec¬ 
tions for heavy nuclei with respect to deuteron for the 
Bjorken scaling parameter £b between about 1.5 and 1.9, 
see, e.g., ref. [28] for a recent review. Systematic anal¬ 
yses of these inclusive experiments and data from ex¬ 
clusive two-nucleon knockout reactions induced by high- 
energy electrons or protons have firmly established that 
the HMT fraction in symmetric nuclear matter (SNM) 
is about = 28% ± 4% and that in PNM is about 

x hmt = 15% ± 0 . 5 % [12-14, 29], 

The Cy|k | 4 shape of the HMT for both SNM and 
PNM is strongly supported by recent findings theoret¬ 
ically and experimentally. The HMT for deuteron from 
variational many-body calculations using several modern 
nuclear forces decrease as |k | -4 within about 10 % and 
in quantitative agreement with that from analyzing the 
d(e, e'p) cross section in directions where final state in¬ 
teraction suffered by the knocked-out proton is small [ 12 ]. 
The extracted magnitude Csnm = Cq of the HMT in 
SNM at po is Co « 0.15 ± 0.03 [12] (properly rescaled 
considering the factor of 2 difference in the adopted nor¬ 
malizations of rik here and that in refs. [12, 29]). Rather 
remarkably, a very recent evaluation of medium-energy 
photonuclear absorption cross sections has also presented 
clear and independent evidence for the C/|k | 4 behavior of 
the HMT and extracted a value of Co « 0.172±0.007 [10] 
for SNM at po in very good agreement with that found 
in ref. [12]. In the following, we use Co ~ 0.161 ± 0.015 
from taking the average of the above two constraints. 
With this Co and the value of given earlier, the 

HMT cutoff parameter in SNM is determined to be 
0o = (1 - x snm/3Co) _1 = 2.38 ± 0.56. 

Very interestingly, the 1/|k | 4 behavior of the HMT nu¬ 
cleons is identical to that in two-component (spin-up and 
-down) cold fermionic atoms first predicted by Tan [30] 
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and then quickly verified experimentally [31]. Tan’s gen¬ 
eral prediction is for all two-component fermion systems 
having an s-wave contact interaction with a scattering 
length a much larger than the inter-particle distance d 
which has to be much longer than the interaction range 
r e . At the unitary limit when |fcFa| —» oo, Tan’s predic¬ 
tion is universal for all fermion systems. Since the HMT 
in nuclei and SNM is known to be dominated by the ten¬ 
sor force induced neutron-proton pairs with the a ~ 5.4 
fm and d « 1.8 fm at p 0 , as noted in refs. [10, 12], Tan’s 
stringent conditions for unitary fermions is obviously not 
satisfied in normal nuclei and SNM. The observed identi¬ 
cal l/|k| 4 behavior of the HMT in nuclei and cold atoms 
may have some deeper physical reasons deserving fur¬ 
ther investigations. Indeed, a very recent study on the 
A(e,e'p) and A(e,e'pp) scattering has shown that the ma¬ 
jority of the short range correlation (SRC)-susceptible 
n-p pairs are in the 3 Si state [32]. On the other hand, 
because of the unnaturally large neutron-neutron scatter¬ 
ing length a nn ( 1 So) = —18.8 fm, it is known that PNM 
is closer to the unitary limit [33]. The EoS of PNM can 
thus be expanded as [34] 


as shown in the inset, the Eq. ( 6 ) is consistent with 
the prediction by the effective field theory [33]. In the 
range of 0.01 fm -3 to about 0 . 02 fm -3 , it has some de¬ 
viations from predictions in ref. [33] but agrees very well 
with the NLO lattice simulations [37]. At higher densities 
up to about po, it overlaps largely with predictions by the 
chiral perturbation theories [38] and the quantum Monte 
Carlo simulations [39, 40]. In addition, recent studies on 
the spin-polarized neutron matter within the chiral effec¬ 
tive field theory including two-, three-, and four-neutron 
interactions indicate that properties of PNM is similar 
to the unitary Fermi gas at least upto po far beyond the 
scattering-length regime of p < po/100 [41]. Overall, the 
above comparison and studies clearly justify the use of 
Eq. ( 6 ) to calculate the PNM EoS up to about p 0 . 

Both the HMT and EoS can be experimentally 
measured independently and calculated simultaneously 
within the same model. Tan has proven in great detail 
that the two are directly related by the so-called adiabatic 
sweep theorem [30]. It is valid for any two-component 
Fermi systems under the same conditions as the Eq. 6 
near the unitary limit. For PNM, it can be written as 


-Epnm(p) 


3 (fcg NM ) 2 r c_ 

5 2 M _ fcp NM a nn 3(fc£ NM a nn ) 2 


( 6 ) 

where fcp NM = 2 1 ' 3 fcF is the transition momentum in 
PNM, £ ss 0.4 ± 0.1 is the Bertsch parameter [35], ( ss 
v k, 1 are two universal constants [36]. 



FIG. 2: (Color Online) The EoS of PNM obtained from 
Eq. (6) (dashed red band) and that from next-leading-order 
(NLO) lattice calculation [37] (blue solid points), chiral per¬ 
turbative theories [38] (green band), quantum Monte Carlo 
simulations (QMC) [39, 40] (magenta band and purple stars), 
and effective field theory [33]. 


Shown in Fig. 2 is a comparison of the EoS of PNM ob¬ 
tained from Eq. ( 6 ) (dashed red band) with several state- 
of-the-art calculations using modern microscopic many- 
body theories. At densities less than about 0.01 fm~ 3 , 


c n PNM • (fcr M ) 4 = -47 tM • . (7) 

d(a 3 ) 

While the results shown in Fig. 2 justify the use of Eq. 6 
for the EoS of PNM up to about po , indeed, to our best 
knowledge there is currently no proof that the Eq. 7 is 
also valid in the same density range as the Eq. 6 . Thus, 
it would be very interesting to examine the validity range 
of Eq. 7 using the same models as those used to calculate 
the EoS. In this work, we assume that the Eqs. 6 and 
7 are both valid in the same density range. Then, the 
strength of the HMT in PNM can be readily obtained as 

» 2C/5tt + 4zz/(37rA:p NM a nn ( 1 S 0 )) « 0.12. (8) 

Noticing that C PNM = Co(l + Ci), we can then infer that 
C\ = —0.25 ± 0.07 with the Co given earlier. Next, after 
inserting the values of ajp^M an d C PNM into Eq. (2), the 
high momentum cutoff parameter for PNM is determined 
to be 0 PNM = M 1 + 0 i) = (1 - = 

1.04 ± 0.02. It is not surprising that the </> PNM is very 
close to unity since only about 1.5% neutrons are in the 
HMT in PNM. Subsequently, using the cf> o determined 
earlier, we get cj> i = —0.56 ± 0.10. 

The two parameters 0o and /?i in fdj = /? 0 (1 + /3it/< 5) 
depend on the function J(|k|//cp) which is still model de¬ 
pendent. To minimize the model assumptions and eval¬ 
uate the dominating terms in the kinetic EoS, in the fol¬ 
lowing we shall first use a momentum-independent deple¬ 
tion of the Fermi sea as in most studies in the literature. 
The HMT parameters Cj and 4>j evaluated above remain 
the same. Then, we examine the maximum correction to 
each term in the kinetic EoS by using the largest values 
of [3 o and /3i allowed and a typical function 7(|k|/fcp). 
Not surprisingly, the corrections are all small. 
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IV. ISOSPIN DEPENDENCE OF KINETIC EOS 
OF ANM 


The kinetic EoS can be expanded in S as 

E kin (p, 5) = E^(p) + E^ m (p) 6 2 + E^ mA (p) 6 4 + 0(S 6 ). 

(9) 

The coefficients evaluated from Eq. (4) using the ? 7 j/(p, <5) 
in Eq. (1) with /3j = 0 are 


Er(p) =jE F (p) 


CW =3^ f (p) 


1 + Cq ( 50O + t-' 


( 10 ) 


1 + Co (1 + 3Ci) ( 50o + —- 


+ 3Co</>i ( 1 + —Ci ) ( 50o — 


00 


+ 


27 C o 0‘ 


50o 


(ID 


E^ mA (p) =^E F (p) 


1 + C 0 (l — 3Ci) ( 50o + —-I 


+ 3Co</> 1 (9C 1 — 1) ( 50o — 


+ 


00 

810001(901 - 9Ci<^i - 15^ + 15Ci + 5) 
50o 


■ ( 12 ) 


In the FFG where there is no HMT, 0o = 1, 0i = 0 and 
thus 50 o+3/0o—8 = 0 , the above expressions reduce nat¬ 
urally to the well known results of Eo(p) = 3E F (p)/5 , 
E^(p) = E F (p)/ 3, and E^ A (p) / E^ m (p) = 1/27 
where E F (p) — kp/2M is the Fermi energy. 

For the interacting nucleons in ANM with the mo¬ 
mentum distribution and its parameters given earlier, we 
found that Eg in (p 0 ) = 40.45 ± 8.15 MeV, £ k ^(p 0 ) = 
-13.90 ± 11.54 MeV and E^l lA (p () ) = 7.19 ± 2.52 MeV, 
respectively. Compared to the corresponding values for 
the FFG, it is seen that the isospin-dependent HMT in¬ 
creases significantly the average kinetic energy £n(po) 
of SNM but decreases the kinetic symmetry energy 
Esym(Po) °f ANM to a negative value qualitatively con¬ 
sistent with findings of several recent studies of the ki¬ 
netic EoS considering short-range nucleon-nucleon corre¬ 
lations using both phenomenological models and micro¬ 
scopic many-body theories [42-47]. However, it was com¬ 
pletely unknown before if the empirical isospin parabolic 
law is still valid for the kinetic EoS of ANM when the 
isospin-dependent HMTs are considered. Very surpris¬ 
ingly and interestingly, our calculations here show clearly 
that it is broken seriously. More quantitatively, the ratio 
|A ,k /m i 4 (po)/A’^m(po)| is about 52% ± 26% that is much 
larger than the FFG value of 3.7%. We also found that 
the large quartic term is mainly due to the isospin de¬ 
pendence of the HMT cutoff described by the 0i param¬ 
eter. For example, by artificially setting 0i = 0, we ob¬ 
tain E^(p 0 ) = 14.68 ± 2.80 MeV and E^ A (p 0 ) = 
1.12 ± 0.27 MeV which are all close to their FFG values. 

Considering short-range nucleon-nucleon correlations 
but assuming that the isospin parabolic approximation 


is still valid, some previous studies have evaluated the 
kinetic symmetry energy by taking the difference 

between the kinetic energies of PNM and SNM, i.e., sub¬ 
tracting the Ep/J M by Eq 111 . This actually approximately 
equals to E k /^(p 0 ) + E^^ ^po) = -6.71 ± 9.11 MeV in 
our current work. This value is consistent quantitatively 
with the E^(po) found in ref. [29] using the parabolic 
approximation. 


V. CORRECTIONS DUE TO THE 
MOMENTUM-DEPENDENT DEPLETION OF 
THE FERMI SEA 


To estimate corrections due to the momentum depen¬ 
dence of the depletion very close to the Fermi surface, i.e., 
a finite /3j, we consider a widely used single-nucleon mo¬ 
mentum distribution parameterized in ref. [25] based on 
calculations using many-body theories. For |k| < 2 fm _ , 
it goes like ~ e _a l k l with a ~ 0.12frn 2 . At po since 
ak F ~ 0 . 21 , e _Q l k l 2 ~ 1 — a|k| 2 + 0(|k| 4 ) is a good 
approximation in the range of 0 < |k| < kp. Thus, 
we adopt a quadratic function J(|k|/fcp) = (|k|/fcp) 2 . 
The constants in the parameterization of ref. [25] are ab¬ 
sorbed into our parameters A j and /3j. Then Eq. (3) 
gives us A,/ = 1 — 3/3j/5 — 3Cj (1 — l/0j). Specifi¬ 
cally, we have 0o = (5/3) [1 — A 0 — 3Co(l — 0Q 1 )] = 
(5/3) [1 — A 0 — £g]^] for SNM. Then using the predicted 
value of Ao ~ 0.88 + 0.03 [9, 22, 23] and the experimental 
value of « 0.28 ± 0.04, the value of /3q is esti¬ 

mated to be about —0.27 ±0.08. Similarly, the condition 
Pj = A)(l + P 1 T 3 S) < 0, i.e., is a decreasing function 
of momentum towards k F , indicates that |/3i| < 1. 
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FIG. 3: (Color Online) Corrections to the E^ n (po) and 
Esym,4,(P°) as functions of /3i with /3o = —0.35. 


First of all, a finite value of (3j is expected to affect 
the “renormalization function” Zp. For SNM, we have 
Z% = l + 2/3 0 /5-C' 0 -a;f/) I J = 0.45 + 0.07 (0.56 + 0.04) in 
the presence (absence) of /3 q- For ANM, however, the Zp 
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depends on the less constrained value of /? 4 . It is worth 
noting that the latter also determines the neutron-proton 
effective E-mass splitting which has significant effects on 
isovector observables in heavy-ion collisions [48], and a 
study is underway to further constrain the value of /3i 
using data from heavy-ion reactions. 

Contributions from a finite /3j to the first three terms 
of the kinetic EoS are 


O) w + 


1 


<^sym,4 (P) = ^ E f{Po) 


4fl)(l ~ 3/30 
35 


(13) 

(14) 

(15) 


With the largest magnitude of /3o = —0.35, we examine 
in Fig. 3 the corrections to the E^^po) and ^ym j4 (po) 
as functions of f3\ in its full range allowed. In this case 
the maximum effects of the finite f3j are revealed. It is 
seen that the correction on the E sym, 4 ,(Po) negligible 
while the correction on the Eg^^po) is less than 2MeV. 
Considering the corrections due to the finite /3o and f3\ 
and their uncertainties, we finally obtain EQ ln (po) = 
39.77 ± 8.13 MeV, E^ (p 0 ) = -14.28 ± 11.59 MeV and 
-®sym, 4 (po) = 7.18 ± 2.52 MeV, respectively. We notice 
here that the 6 6 term was also consistently evaluated and 
was found to be negligibly small at po- 


VI. SUMMARY AND DISCUSSIONS 

In summary, using an isospin-dependent single-nucleon 
momentum distribution including a high (low) momen¬ 
tum tail (depletion) with its shape parameters con¬ 
strained by the latest results of several relevant experi¬ 
ments and the state-of-the-art predictions of modern mi¬ 
croscopic many-body theories, we found for the first time 
that the kinetic EoS of interacting nucleons in ANM is 
not parabolic in isospin asymmetry. It has a significant 
quartic term of 7.18 ± 2.52 MeV while its quadratic term 
is —14.28 ± 11.60 MeV at saturation density of nuclear 
matter. 

To this end, it is necessary to point out the limitations 
of our approach and a few physical implications of our 
findings. Since we fixed the parameters of the nucleon 
momentum distribution (Eq. (1)) by using experimental 
data and/or model calculations at the saturation density, 
the possible density dependence of these parameters is 
not explored in this work. The density dependence of the 
various terms in the kinetic EoS is thus only due to that 
of the Fermi energy as shown in Eqs.(10)-(12). In this 
limiting case, the slope of the kinetic symmetry energy, 
i.e., L kin = 3p 0 &E^£>)/dp| p=Po = -27.81 ± 23.08 MeV 
while that of the FFG is about 25.04 MeV. 


The SRC-reduced kinetic symmetry energy with re¬ 
spect to the FFG prediction has been found to affect 
significantly not only our understanding about the origin 
of the symmetry energy but also several isovector ob¬ 
servables, such as the free neutron/proton and n~/ tt + 
ratios in heavy-ion collisions [29, 49, 50]. However, to 
our best knowledge, an investigation on possible effects 
of a large isospin quartic term on heavy-ion collisions has 
never been done while its effects on properties of neutron 
stars have been studied extensively [3]. Of course, effects 
of the quartic and quadratic terms should be studied to¬ 
gether within the same approach. To extract from nu¬ 
clear reactions and neutron stars information about the 
EoS of neutron-rich matter, people often parameterize 
the EoS as a sum of the kinetic energy of a FFG and 
a potential energy involving unknown parameters upto 
the isospin-quadratic term only. Our findings in this 
work indicate that it is important to include the isospin- 
quartic term in both the kinetic and potential parts of 
the EoS. Moreover, to accurately extract the completely 
unknown isospin-quartic term E^°^ 4 (p)d 4 in the poten¬ 
tial EoS it is important to use the kinetic EoS of quasi¬ 
particles with reduced kinetic symmetry energy and an 
enhanced quartic term due to the isospin-dependence 
of the HMT. Most relevant to the isovector observables 
in heavy-ion collisions, such as the neutron-proton ra¬ 
tio and differential flow, is the nucleon isovector poten¬ 
tial. Besides the so-called Lane potential ±2pEF°^(p)S 
where the (p) is the potential part of the symme¬ 
try energy and the ± sign is for neutrons/protons, the 
-Esym 4 (p)<5 4 term contributes an additional isovector po¬ 
tential ±4pE^°^ 4 (p)S 3 . In neutron-rich systems besides 
neutron stars, such as nuclear reactions induced by rare 
isotopes and peripheral collisions between two heavy nu¬ 
clei having thick neutron-skins, the latter may play a sig¬ 
nificant role in understanding the isovector observables 
or extracting the sizes of neutron-skins of the nuclei in¬ 
volved. We plan to study effects of the isospin-quartic 
term in the EoS in heavy-ion collisions using the isospin- 
dependent transport model [4] in the near future. 
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